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The Emilia seismic sequence (Northern Italy) started on May 2012 and caused 17 casualties, severe damage 
to dwellings and forced the closure of several factories. The total number of events recorded in one month 
was about 2100, with local magnitude ranging between 1.0 and 5.9. We investigate potential mechanisms 
(static and dynamic triggering) that may describe the evolution of the sequence. We consider rupture 
directivity in the dynamic strain field and observe that, for each main earthquake, its aftershocks and the 
subsequent large event occurred in an area characterized by higher dynamic strains and corresponding to 
the dominant rupture direction. We find that static stress redistribution alone is not capable of explaining 
the locations of subsequent events. We conclude that dynamic triggering played a significant role in driving 
the sequence. This triggering was also associated with a variation in permeability and a pore pressure 
increase in an area characterized by a massive presence of fluids. 

The Emilia region, Northern Apennines, Italy, was struck by a seismic sequence that started on May 19, 2012, 
at 23:13:27 GMT with a Mw 3.8 earthquake. It produced about 2100 events during the following month, 
affecting an area of about 60 km X 30 km elongated in the EW direction (Figure 1). The sequence began 
with the reactivation of two buried, sub-parallel N100 c E-striking thrust faults with hypocenters located mainly in 
the upper 10 km 1 . The largest events, Mw 5.6 and 5.4, occurred on May, 20 and 29, respectively, and were followed 
by 6 Mw > 4.5 earthquakes over the next 2 weeks. 

This is a strategic area for the Italian economy and a site of intensive petroleum extraction 2 and gas storage 
(http://unmig.sviluppoeconomico.gov.it/unmig/pozzi/completo.asp). For this reasons, understand the physical 
processes involved in triggering of this seismic sequence, with particular attention on the role of fluids, and try to 
improve the risk assessment is particularly important in this part of Italy. In this study we investigate the static 
effect of the Coulomb stress redistribution and some dynamic effects of passing seismic waves. 



Results 

We use 22 hypocentral locations, moment magnitudes and focal mechanism solutions 3 to estimate source 
dimensions 4-5 and cumulative changes in the static stress field (ACFS). In addition, we use peak-ground velocities 
(PGVs) of the 8 largest earthquakes (Table 1) to estimate: i) rupture directivity, and ii) the peak-dynamic strain 
field. Combining i) and ii) results in an original representation of the dynamic strain field, whose amplitude and 
azimuthal distribution is modified by source directivity. The PGVs used in this study are those also used to 
produce ShakeMaps in Italy 6 . The waveforms are recorded by the Italian permanent seismic network (http:/ /iside. 
rm.ingv.it) and Rete Accelerometrica Nazionale (http://www.protezionecivile.gov.it/jcms/it/ran.wp). The com- 
bination of the two datasets provides data with a reliable azimuthal coverage of stations. According to the 
automatic procedures implemented at INGV 6 , seismograms are corrected by the instrumental response, and 
processed applying a de-trending and a band-pass filtering in the range 0.01-30 Hz. At each station, the PGV 
corresponds to the largest value between the two horizontal components of the recorded velocity. The number of 
PGV data available for each earthquake is reported in Figure 2 where only those seismic stations located at 
epicentral distances less than 150 km are included. Moreover, we discarded peak velocity values differing more 
than 2c> from the predictions provided by the Ground-Motion Prediction Equation (GMPE) 7 , where cr is the 
standard error of the GMPE. This resulted in discarding about 18% on average, of available data. 
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Figure 1 | Location map of the 2012 Emilia earthquake sequence. Geographic location of the 2012 Emilia seismic sequence (Italian Seismological 
Instrumental and parametric Data-basE, http://iside.rm.ingv.it). Orange dots, with corresponding focal mechanisms 3 , represent the epicenters of the 
main events analysed in this study and listed in Table 1 . Red dots represent the location of events 2 and 7, whose focal mechanisms have been estimated by 
average of the known focal mechanisms. For each event the sequential number and magnitude are also specified. Open stars indicate location of events 
with M > 3.0, while gray squares represent all the other earthquakes. The sequence spans a time interval of about one month and covers an area of about 
60 km X 30 km. The figure was generated by using the Generic Mapping Tools (http://gmt.soest.hawaii.edu/) 35 . The topography is based on the 
ETOPOl 1 Arc-Minute Global Relief Model, http://www.ngdc.noaa.gov/mgg/global/global.html 3 ''. 



The adopted GMPE is obtained from the analysis of the Italian 
earthquakes and is valid for magnitude ranging between 4.0 and 6.9 
and distances up to 200 km (ref. 7). 

From the analysis of the cumulative stress transfer (ACFS) we 
argue that its effect at the locations and on the focal mechanisms of 
the largest subsequent earthquakes does not explain their occur- 
rence. Indeed, for the first five main events, the corresponding 
ACFS values are one order of magnitude smaller than the minimum 
ACFS commonly required to significantly contribute to the trigger- 
ing process alone (i.e., 0.01 MPa) 8,9 (supplementary Figure SI and 
also Figure 3a). This statement holds regardless of whether we 
consider the preferential nodal plane (inferred from geological 
information, i.e. E-W striking, S-dipping 110 ) or the auxiliary one. 
Additionally, since static stress redistribution is commonly charac- 
terized by symmetrical lobes of positive and negative values around 
the seismic source, it makes difficult to identify any preferred dir- 
ection for subsequent event locations. Moreover, the symmetry of 



the static stress fields also differs from the asymmetries in the after- 
shock patterns. 

Summarizing, although static stress changes may affect the evolu- 
tion of this sequence, they not appear to be significant contributors to 
the triggering process. In the following of our study, we try to over- 
come the aforementioned limitation of the static Coulomb model by 
combining the source rupture directivity and the dynamic strain 
field. 

Previous studies 1119 also have noted the influence of the rupture 
direction in the dynamic stress field for large and small-to-moderate 
earthquakes. 

We analysed the major events that occurred during the sequence 
(Table 1 and Figure 1) by considering only those earthquakes for 
which focal mechanisms and/or aftershock distributions are avail- 
able. We identified the horizontal rupture directions 4> and estimated 
the Mach numbers a (i.e., the ratio between rupture and shear-wave 
velocities) through the inversion of PGVs 20 ' 21 . We show in Figure 2 



Table 1 | Source parameters of the main earthquakes of the Emilia seismic sequence 

Hypocentral location 



Event 


Date 


Time 


Lat(N) 


Long(E) 


Depth(km) 


Mw 


4> 


5 


X 


<t 


D 


4 


>s 


1 


05/20/2012 


02:03:52 


44.889 


1 1.228 


6.0 


5.6 


105/285 


45/45 


90/90 


102 


± 12 


286 


± 12 


2* 


05/20/2012 


02:07:31 


44.863 


1 1 .370 


5.0 


5.1 


77/252 


37/53 


94/87 


298 


± 35 


98 


± 49 


3 


05/20/2012 


03:02:50 


44.860 


1 1 .095 


13.0 


4.7 


132/245 


50/65 


147/45 


80 


± 24 






4 


05/20/2012 


13:18:02 


44.831 


1 1 .490 


7.0 


4.8 


1 20/275 


38/55 


1 1 0/75 


292 


± 26 


110 


± 21 


5 


05/29/2012 


07:00:03 


44.851 


1 1.086 


5.0 


5.4 


97/270 


45/45 


95/85 


304 


± 22 


116 


± 20 


6 


05/29/2012 


10:55:57 


44.888 


1 1 .008 


6.0 


5.1 


96/282 


40/50 


85/94 


278 


± 73 






7* 


05/29/2012 


1 1:00:25 


44.8790 


1 0.9470 


5.4 


5.2 


77/252 


37/53 


94/87 


106 


±41 






8 


06/03/2012 


19:20:43 


44.899 


10.943 


10.0 


4.6 


1 17/265 


29/65 


1 1 9/75 


76 


± 24 


208 


± 38 



Hypocentral location, moment magnitude, and fault plane solution (strike, (j>, dip, 5, and rake, X, for the principal and auxiliary planes) of the main earthquakes of the Emilia seismic sequence 3 . For the events 
marked by an asterisk, locations are taken from the INGV catalog (http://iside.rm.ingv.it) and focal mechanisms are an average of the all 22 known focal mechanisms 3 . The last two columns list the 
dominant (4>d) and secondary ((j> s ) rupture directions, respectively, estimated from the PGVs inversion. The uncertainties on 4>D an d 4>s have been evaluated measuring the Half Width at Half Maximum of the 
best-model Pdfs shown in Figure S2. 
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Figure 2 | Maps of dominant rupture direction. Overview of inferences obtained from the inversion of the peak-ground velocities for the main 
triggering events (black stars) considered in this study. In each panel, the dominant rupture directions (red arrows) and aftershock distributions 
(gray stars) up-to the next main event in the sequence (green stars), are shown. The orange arrow indicates the secondary rupture direction and its length 
depends on the ratio between the Pdf at the secondary maximum and the value at the absolute maximum (supplementary Figure S2) . The uncertainties of 
the rupture directions are given in Table 1. The blue arrows indicate strike directions of the principal and auxiliary fault planes deduced from the 
corresponding focal mechanism 3 . The black triangles identify the stations at which PGVs data were available. In each panel the origin time and magnitude 
of the triggering event are specified, as well as the best value of the Mach number a and the available number of PGV data used. The figure was generated by 
using the Generic Mapping Tools (http://gmt.soest.hawaii.edu/) 35 . 



locations of triggering events, their dominant rupture directions and 
aftershocks distribution occurred up-to the next event in the 
sequence. From a visual inspection, a correlation is evident: we 
observe that for each main earthquake its aftershocks and the sub- 
sequent main event occur in the areas oriented as the source dom- 
inant rupture direction. 

The probability density functions (Pdfs) of c|) (supplementary 
Figure S2) are analysed to assess the reliability of the retrieved rup- 
ture direction. The uncertainties, evaluated as the Half Width at Half 
Maximum of the Pdfs corresponding to the best model, are of the 



order of 10-40 degrees (Table 1). In some cases (e.g., events 4, 5 and 
8) the 4>-Pdfs present a secondary mode that indicates a deviation 
from a purely unilateral rupture. Note that we consider the secondary 
maximum only if it differs more than 90 degrees from the principal 
one. For the aforementioned events the Pdfs on the percent unilateral 
rupture parameter e (see Method section for its formal definition) also 
suggest a quasi-bilateral rupture propagation on the fault plane. It is 
worth noting that we find different best-fit Mach number a that either 
0.6 or 0.9. This result may be attributed to a variation in the rupture 
velocity rather than a variation in the shear-wave velocity. Indeed, 
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Figure 3 | Static and dynamic stress transfer and permeability change. Panel (a): local cumulative ACFS at hypocentral locations and on preferential 
focal mechanisms of target main events. Cumulative ACFS is estimated considering the contribution of all the earthquakes occurred before the target 
event. The dataset used for Coulomb stress estimation contains 22 events for which locations and focal mechanisms are available 3 . For the events 2 and 7 in 
Table 1 the fault plane solutions are not available and a mean focal mechanism is then assumed. The dashed line represents the commonly accepted 
triggering threshold for static ACFS 8,9 . Below this threshold, stress changes are considered too small to significantly contribute to the triggering process. 
Panel (b): local dynamic stress change obtained from the peak dynamic strain (assuming a rigidity value of 30 GPa) induced by each considered event at 
the hypocenter of the next main earthquake in the sequence. The local dynamic stress is estimated both considering the directivity function Cd, (squares) 
and ignoring it (open circles). Panel (c): permeability change AK induced by each considered event at the hypocenter of the next main earthquake in the 
sequence. The two horizontal dashed lines define the seismogenic permeability range (5 X 10~ 16 - 5 X 10~ 14 m 2 ) 27,28 . AK changes are estimated 
through PGV values: i) modified by the source directivity function Cd (squares); and ii) not modified by Cd (open circles). The gray triangles represent 
AK obtained by using the upper and lower R bounds (see equation 2). The relative hypocentral distance (in km) between each triggering and target events, 
is provided in the bottom of the panel. In the three panels, the numbering on x-axis refers to the seismic events as listed in Table 1. 



rupture velocity variations produce a different ground-motion field: 
the faster the rupture, the larger the ground-motion amplitudes 22 24 . 

PGVs can be also used to estimate the peak-dynamic strain field 18 . 
We use the directivity function (see Method section) to account 
for modifications of PGVs due to source directivity while computing 
the dynamic strain field. 



Figure 4 shows that the major events of the Emilia seismic 
sequence occurred in areas where the peak-dynamic strain values 
are in the range of the dynamic triggering 12 (from a few microstrain 
up to tens of microstrain). Moreover, the corresponding local dyna- 
mic stress changes (—0.1-1 MPa, assuming a rigidity of 30 GPa) are 
about one order of magnitude higher than the cumulative static 
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Figure 4 | Geographical distribution of peak-dynamic strain. Peak-dynamic strain field (microstrain) obtained by using the peak-ground velocity 
as a proxy (e ~ PGV/Vs). In each panel, the black and green stars indicate the triggering and target event, respectively. The red and orange arrows represent 
the dominant and the secondary rupture directions, respectively. Black rectangles represent the surface fault projections estimated through PGVs 
inversion (in some cases dimensions are very small and then dimly visible). 



stress changes (Figure 3b and 3a). This might suggest that the 
dynamic effect played a significant role in the evolution of this seismic 
sequence. 

Furthermore, we analyse the variation of the permeability, Ak, in 
the fractured area. Permeability variations are also associated with 
the passage of seismic waves, which indeed increase pore pressure 
and fluid-flow rates in a fluid-saturated medium 25 . The presence of 
fluids in the Emilia region has been suggested in a tomographic study 
that found a high Vp and high Vp/Vs ratio (about 1.814) for the area, 
which indicates a fluid-saturated and partially-fractured medium 26 . 
Our idea is that source directivity can also enhance changes in per- 
meability. Such variations of permeability might encourage the frac- 
ture in a preferred direction as a consequence of a local increase of the 
pore pressure and fluid flow rates, resulting in a reduction of the 
effective normal stress. For each main event we evaluate the change 
in permeability at the location of the next earthquake in the sequence. 
For all the considered events, the induced Ak values are in the range 
of the seismogenic permeability 27,28 (from 5 X 10~ 16 m 2 to 5 X 
10~ 14 m 2 ) (Figure 3c) and can then be associated with the evolution 
of the seismic sequence. 

Discussion 

This study indicates that the dynamic triggering caused by passing 
seismic waves and enhanced by source directivity might be the 



primary factor to explain the evolution of the 2012 Emilia seismic 
sequence. In fact, we observed a correlation between the locations of 
aftershocks and subsequent main events with: i) the peak dynamic 
strain fields; ii) the local change of the permeability. 

We also observe that static Coulomb stress changes did not play a 
significant role, remaining under the threshold value for static trig- 
gering. 

The correlation between variations of permeability and location of 
events suggests that the presence of fluids might contribute to the 
evolution of the sequence by altering the effective normal stress in 
specific areas. 

We think that this issue deserves future research to be used to 
improve the forecasting ability of models based on stress change. 

Methods 

We estimate the surface fault projection and the dominant horizontal rupture dir- 
ection of the analysed earthquakes through an inversion technique of PGVs. This 
technique minimizes the difference between observed and predicted PGV values 
using a grid-searching scheme in a Bayesian framework 20 . Peak predictions are 
obtained from a GMPE by using the Rf B distance (i.e., the minimum distance of a site 
from the surface fault projection) metric 29 . A first order correction for site effects is 
obtained through a method 30 that uses the value of the average shear-wave velocity 
over the uppermost 30 m (VS30) and provides with multiplicative coefficients. The 
VS30 values have been extracted from the database compiled for implementing 
ShakeMap in Italy 6 . Predictions including site effect are then modified to account for 
the source directivity through the directivity function given by 
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where .9 is the angle between the ray leaving the source and the direction of rupture 
propagation <J) 31 , and a is the Mach number. The percent unilateral rupture e, is 
defined as {21' — L)/L, where L is the fault length and V the prevalent rupture 
direction 17 . It allows for the discrimination between unilateral (e — 1) and bilateral (e 
— 0) ruptures. 

Details on the inversion technique are provided in the paper in which the method 
has been originally presented 20 . Here, we just note that the best solution for surface 
fault projection, and for <J) and e, corresponds to the largest value of the a-posteriori 
Pdfs resulting from the adopted Bayesian approach. The associated standard error are 
obtained by analysing the marginal Pdfs. 

As we have analysed a sequence of earthquakes characterized by different (mod- 
erate) magnitude, we expect that a might not be the same for all the events. Thus, for 
each earthquake we inverted PGV values by exploring a in the range 0.2-0.9 with a 0. 1 
increment. The smallest residual value allows us to select the best model. 

The estimation of the peak dynamic strain (e) field is based on an empirical 
approach that uses PGVs as a proxy: e ~ PGV/Vs, where Vs is the shear-wave 
velocity 18 . PGVs predicted by using the GMPE are here modified by the Cd coefficient 
and divided by a factor of 2 in order to correct for the free- surface amplification 32 . 

We estimate variations in the permeability factor through the relationship 25 : 



Ak = R 



PGVC d 
Vs ' 



(2) 



The coefficient R accounts for the strain effect on the permeability response and 
ranges between 3.0 X 10" 10 m 2 and 8.4 X 10" 10 m 2 (ref. 25) with a central value of 5.7 
X 10~ 10 m 2 , which has been used in our study. As for the strain field computation, the 
PGVs account for the directivity function Cd and free-surface amplification. Both for 
peak dynamic strain field and for permeability factor computation, the Vs values at 
each earthquake hypocenter have been estimated from a tomographic velocity model 
valid for the Emilia region 26 . 

Coulomb static stress variations, ACFS, are computed by using the solutions for a 
homogeneous elastic half-space 33 . The Coulomb model is commonly defined as ACFS 
= At + jj.' Aa, where At is the shear stress change on the failure plane (positive in the 
slip direction), Aa is the effective normal stress change (positive for extension) and u' 
is the apparent friction coefficient 34 . We estimated dimensions and mean slip of 
extended square fault patches by using scaling relationships 4,5 and assuming a con- 
stant stress drop of 2.4 MPa on the faults, which is approximately the average value 
retrieved for the study area from the analysis of strong-motion records 35 . 

In the present study we used a friction coefficient u = 0.8 (common to almost all 
the type of rocks), a Skempton ratio B — 0.5, and a rigidity value of 30 GPa. 
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